home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Ham Radio 2000
/
Ham Radio 2000.iso
/
ham2000
/
satellit
/
vstsrc
/
compute.c
< prev
next >
Wrap
C/C++ Source or Header
|
1995-01-23
|
29KB
|
871 lines
/*
* %W% %E% %U% [EXTREL_1.2]
*
* VersaTrack orbit calculations are based on those that in Dr. Manfred
* Bester's sattrack program (the Unix(tm) versions 1 & 2).
*
* The data from which the maps where generated come from "xsat", an
* X-Windows program by David A. Curry (N9MSW).
*
* Site coordinates come from various sources, including a couple of
* World Almanacs, and also from both of the programs mentioned above.
*
* The following are authors' applicable copyright notices:
*
*
* Copyright (c) 1992, 1993, 1994 Manfred Bester. All Rights Reserved.
*
* Permission to use, copy, modify, and distribute this software and its
* documentation for educational, research and non-profit purposes, without
* fee, and without a written agreement is hereby granted, provided that the
* above copyright notice and the following three paragraphs appear in all
* copies.
*
* Permission to incorporate this software into commercial products may be
* obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,
* Berkeley, CA 94709, USA.
*
* IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
* SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF
* THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED
* OF THE POSSIBILITY OF SUCH DAMAGE.
*
* THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
* PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
* BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT,
* UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
*
*
* [xsat for X-windows] Copyright 1992 by David A. Curry
*
* Permission to use, copy, modify, distribute, and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice appear in all copies and that both that copyright
* notice and this permission notice appear in supporting documentation. The
* author makes no representations about the suitability of this software for
* any purpose. It is provided "as is" without express or implied warranty.
*
* David A. Curry, N9MSW
* Purdue University
* Engineering Computer Network
* 1285 Electrical Engineering Building
* West Lafayette, IN 47907
* davy@ecn.purdue.edu
*
* VersaTrack Copyright 1993, 1994 Siamack Navabpour. All Rights Reserved.
*
* Permission is hereby granted to copy, modify and distribute VersaTrack
* in whole, or in part, for educational, non-profit and non-commercial use
* only, free of charge or obligation, and without agreement, provided that
* all copyrights and restrictions noted herein are observed and followed, and
* additionally, that this and all other copyright notices listed herein
* appear unaltered in all copies and in all derived work.
*
* This notice shall not in any way void or supersede any of the other authors
* rights or privileges.
*
* VersaTrack IS PRESENTED FREE AND "AS IS", WITHOUT ANY WARRANTY OR SUPPORT.
* YOU USE IT AT YOUR OWN RISK. The author(s) shall not be liable for any
* direct, indirect, incidental, or consequential damage, loss of profits or
* other tangible or intangible losses or benefits, arising out of or related
* to its use. VersaTrack carries no warranty, explicit or implied, including
* but not limited to those of merchantablity and fitness for a particular
* purpose.
*
* Siamack Navabpour, 12342 Hunter's Chase Dr. Apt. 2114, Austin, TX 78729.
* sia@bga.com or sia@realtime.com.
*/
/* THE FOLLOWING ADAPTED FROM Sattrack (the Unix version by Dr. M. Bester) */
#include <windows.h>
#include <math.h>
#include "vstdefs.h"
#include "vsttype.h"
#include "constant.h"
#include "vstextrn.h"
double
absol( x, y, z)
double x, y, z;
{
return sqrt(x*x + y*y + z*z);
}
double
reduce(value, rangeMin, rangeMax)
double value, rangeMin, rangeMax;
{
double range, rangeFrac, fullRanges, retval;
range = rangeMax - rangeMin;
rangeFrac = (rangeMax - value) / range;
modf(rangeFrac,&fullRanges);
retval = value + fullRanges * range;
if (retval > rangeMax)
retval -= range;
return retval;
}
static double
getkep(ma, eccentricity, keperr)
double ma, eccentricity;
int *keperr;
{
double eccAnomaly; /* eccentric anomaly */
double ta;
register double error;
int count;
count = 5000;
eccAnomaly = ma; /* Initial guess */
*keperr = 0;
do {
error = (eccAnomaly - eccentricity * sin(eccAnomaly) - ma) /
(1.0 - eccentricity * cos(eccAnomaly));
eccAnomaly -= error;
} while (ABS(error) >= EPSILON && --count > 0);
if (count <= 0)
*keperr = 1;
if (ABS(eccAnomaly - PI) < EPSILON)
ta = PI;
else
ta = 2.0 * atan(sqrt((1.0 + eccentricity) /
(1.0 - eccentricity)) * tan(eccAnomaly / 2.0));
return reduce (ta, (double)0.0, TWOPI);
}
/*
* subsatpos: calculates sub-satellite point
*
*
* EARTHRADIUS = equatorial radius of the Earth (m)
* from "Astronomical Almanac", 1984, p. S7
*
* F,FF,FFF = flattening factors for the Earth (Hayford's spheroid)
* from "Astronomical Almanac", 1984, p. S7
* F1,2,3 = coefficients for the latitude correction formula
* from "Explanatory Supplement to the Ephemeris", 1961, p.58
*
* REFF0,1,2,3 = coefficient for the Earth's radius correction formula
* from "Explanatory Supplement to the Ephemeris", 1961, p.58
*
* effEarthRad = effective radius of the Earth at the location of the
* sub-satellite point with respect to sea level (m)
*
* for reference see: a) "Methods of Experimental Physics", Vol. 12C, p.308
* b) "Explanatory Supplement to the Ephemeris",1961, p.58
*/
static void
ssp(tp, satX, satY, satZ, time, latitude, longitude, height, avis)
track_t *tp;
double satX, satY, satZ, time;
float *latitude, *longitude;
double *height, *avis;
{
double x, y, z, r, b2, b4, b6, lgt, effrad;
r = absol(satX, satY, satZ);
x = satX / r;
y = satY / r;
z = satZ / r;
if (z > 1.0) z = 1.0;
if (z < -1.0) z = -1.0;
*latitude = (float) asin(z);
b2 = 2.0 * (*latitude);
b4 = b2 + b2;
b6 = b2 + b4;
*latitude += (float) ((F1 * sin(b2)) + (F2 * sin(b4)) + (F3 * sin(b6)));
#ifdef SIMPLER_MODEL
lgt = TWOPI * ((time - tp->sidday) * SIDSOLAR + tp->sidtime)
- atan2(y,x);
#else
lgt = tp->gastime - atan2(y,x);
#endif
*longitude = (float) reduce(lgt, -PI, PI);
effrad = (REFF0 + REFF1 * cos(b2) + REFF2 * cos(b4) + REFF3 * cos(b6))
* EARTHRADIUS;
*height = r - effrad;
/* XXX FIX - sia: set crash flag if height <= 0.0 */
*avis = acos(effrad/r); /* geocenteric angle of visibility */
}
/*
* precess: calculates precession of the satellite's orbit
*/
static void
precess(semiMajorAxis, eccentricity, inclination, RAANPrec, perigeePrec)
double semiMajorAxis, eccentricity, inclination;
double *RAANPrec, *perigeePrec;
{
double factor;
/*
* See the Satellite Experimenter's Handbook, pp. 11-11,11-12,11-13.
* The 3rd mulitiplicative term in the formula is reduced to zero
* when orbit is near circluar. i.e., eccentricity almost = 0.
* Equation 11.16.
*/
factor = pow(EARTHRADIUS / semiMajorAxis, 3.5)
/ SQR(1.0 - SQR(eccentricity)) * CDR;
*RAANPrec = 9.95 * factor * cos(inclination);
/*
* See Satellite Experimenter's Handbook, p. 11-9. Equation 11.13a
*/
*perigeePrec = 4.97 * factor * ( 5.0 * SQR( cos(inclination) ) - 1.0 );
}
/*
* satpos : calculates satellite position and velocity in the mean
* equatorial coordinate system (RA, Dec)
*/
static int
satpos(epochTime, epochRAAN, epochArgPerigee, semiMajorAxis,
inclination, eccentricity, RAANPrec, perigeePrec, time, trueAnomaly,
curRaan, argPerigee, X, Y, Z, radius, vX, vY, vZ)
double epochTime, epochRAAN, epochArgPerigee;
double semiMajorAxis, inclination, eccentricity;
double RAANPrec, perigeePrec, time, trueAnomaly;
double *curRaan, *argPerigee, *X, *Y, *Z, *radius, *vX, *vY, *vZ;
{
double RAAN;
double Xw, Yw, vXw, vYw; /* in orbital plane */
double tmp;
double Px, Py, Pz, Qx, Qy, Qz; /* Escobal transformation #31 */
double cosArgPerigee, sinArgPerigee;
double cosRAAN, sinRAAN, cosInclination, sinInclination;
/* See Satellite Experimenter's Handbook p. 11-6. */
*radius = semiMajorAxis * (1.0 - SQR(eccentricity))
/ (1.0 + (eccentricity * cos(trueAnomaly)));
if (*radius < EARTHRADIUS) {
#ifdef _DEBUG_
diag("satellite has crashed!\n");
#endif /* _DEBUG_ */
usermsg("Satellite not in orbit! Elementset may be old or incorrect");
return 0;
}
Xw = *radius * cos(trueAnomaly);
Yw = *radius * sin(trueAnomaly);
tmp = sqrt(GMSGP / (semiMajorAxis * (1.0 - SQR(eccentricity))));
vXw = -tmp * sin(trueAnomaly);
vYw = tmp * (cos(trueAnomaly) + eccentricity);
*argPerigee = epochArgPerigee + (time - epochTime) * perigeePrec;
*curRaan = RAAN = epochRAAN - (time - epochTime) * RAANPrec;
cosRAAN = cos(RAAN);
sinRAAN = sin(RAAN);
cosArgPerigee = cos(*argPerigee);
sinArgPerigee = sin(*argPerigee);
cosInclination = cos(inclination);
sinInclination = sin(inclination);
Px = cosArgPerigee * cosRAAN - sinArgPerigee * sinRAAN * cosInclination;
Py = cosArgPerigee * sinRAAN + sinArgPerigee * cosRAAN * cosInclination;
Pz = sinArgPerigee * sinInclination;
Qx = -sinArgPerigee * cosRAAN - cosArgPerigee * sinRAAN * cosInclination;
Qy = -sinArgPerigee * sinRAAN + cosArgPerigee * cosRAAN * cosInclination;
Qz = cosArgPerigee * sinInclination;
*X = Px*Xw + Qx*Yw; /* Escobal transformation #31 */
*Y = Py*Xw + Qy*Yw;
*Z = Pz*Xw + Qz*Yw;
*vX = Px*vXw + Qx*vYw; /* satellite velocity */
*vY = Py*vXw + Qy*vYw;
*vZ = Pz*vXw + Qz*vYw;
return 1;
}
/*
* sitepos: computes the site postion and velocity in the mean
* equatorial coordinate system. The matrix 'siteMat' is
* used by geototopo to convert geocentric coordinates
* to topocentric (observer-centered) coordinates.
*/
static void
sitepos(tp, siteLat, siteLong, siteAlt, siteTime,
siteX, siteY, siteZ, siteVX, siteVY)
track_t *tp;
double siteLat, siteLong, siteAlt, siteTime;
double *siteX, *siteY, *siteZ, *siteVX, *siteVY;
{
double sb2, sb4, sb6, siteRA, cosRA, sinRA, effSiteRad;
double G1, G2, cosLat, sinLat;
sb2 = 2.0 * siteLat;
sb4 = sb2 + sb2;
sb6 = sb2 + sb4;
#ifdef SIMPLER_MODEL
siteLat += F1 * sin(sb2) + F2 * sin(sb4) + F3 * sin(sb6);
#endif
cosLat = cos(siteLat);
sinLat = sin(siteLat);
effSiteRad = (REFF0 + REFF1 * cos(sb2) + REFF2 * cos(sb4) +
REFF3 * cos(sb6)) * EARTHRADIUS + siteAlt;
G1 = effSiteRad * cosLat;
G2 = effSiteRad * sinLat;
#ifdef SIMPLER_MODEL
siteRA = TWOPI * ((siteTime - tp->sidday) * SIDSOLAR + tp->sidtime) - siteLong;
#else
siteRA = tp->lastime;
#endif
cosRA = cos(siteRA);
sinRA = sin(siteRA);
*siteX = G1 * cosRA;
*siteY = G1 * sinRA;
*siteZ = G2;
*siteVX = -SIDRATE * (*siteY);
*siteVY = SIDRATE * (*siteX);
tp->sitemat[0][0] = sinLat * cosRA;
tp->sitemat[0][1] = sinLat * sinRA;
tp->sitemat[0][2] = -cosLat;
tp->sitemat[1][0] = -sinRA;
tp->sitemat[1][1] = cosRA;
tp->sitemat[1][2] = 0.0;
tp->sitemat[2][0] = cosRA * cosLat;
tp->sitemat[2][1] = sinRA * cosLat;
tp->sitemat[2][2] = sinLat;
}
/*
* satange: calculates satellite range
*/
#ifdef RANGE_USEROUTINE
void satrange(siteX,siteY,siteZ,siteVX,siteVY,
satX,satY,satZ,satVX,satVY,satVZ,range,rangeRate)
double siteX, siteY, siteZ, siteVX, siteVY;
double satX, satY, satZ, satVX, satVY, satVZ;
double *range, *rangeRate;
{
double dX, dY, dZ;
dX = satX - siteX;
dY = satY - siteY;
dZ = satZ - siteZ;
*range = absol(dX ,dY dZ);
*rangeRate = ((satVX-siteVX)*dX + (satVY-siteVY)*dY + satVZ*dZ) / *range;
}
#endif /* RANGE_USEROUTINE */
/*
* geototopo: converts from geocentric mean equatorial coordinates to
* topocentric coordinates
*/
static void
geototopo(satX, satY, satZ, siteX, siteY, siteZ, siteMatT, X, Y, Z)
double satX, satY, satZ, siteX, siteY, siteZ;
double *X, *Y, *Z;
double siteMatT[3][3];
{
satX -= siteX;
satY -= siteY;
satZ -= siteZ;
*X = siteMatT[0][0]* satX + siteMatT[0][1]* satY + siteMatT[0][2]* satZ;
*Y = siteMatT[1][0]* satX + siteMatT[1][1]* satY + siteMatT[1][2]* satZ;
*Z = siteMatT[2][0]* satX + siteMatT[2][1]* satY + siteMatT[2][2]* satZ;
}
/*
* ae: calculates azimuth and elevation
*/
static void
ae(satX, satY, satZ, siteX, siteY, siteZ, siteMatA, azimuth, elevation)
double satX, satY, satZ, siteX, siteY, siteZ;
double siteMatA[3][3];
float *azimuth, *elevation;
{
double x, y, z, r;
geototopo(satX, satY, satZ, siteX, siteY, siteZ, siteMatA, &x, &y, &z);
r = absol(x ,y ,z);
x /= r;
y /= r;
z /= r;
if (x == ZERO)
*azimuth = (y > ZERO) ? (float) HALFPI : (float) THREEHALFPI;
else
*azimuth = (float) (PI - atan2(y, x));
if (z > 1.0) z = 1.0; /* protect against out-of-range problems */
if (z < -1.0) z = -1.0;
*elevation = (float) asin(z);
}
/*
* isinsun: calculates the eclipses of the satellite
*
* this function assumes that the Sun is point-like rather than extended
* and that the Earth is perfectly spherical
* XXX sia - DOUBLE CHECK FOLLOWING.
*/
static int
insun(tp, satX, satY, satZ, siteX, siteY, siteZ, eclTime, satElevation)
track_t *tp;
double satX, satY, satZ, siteX, siteY, siteZ, eclTime;
float satElevation;
{
double meanAnomaly, trueAnomaly;
double sunX, sunY, sunZ, sunVX, sunVY, sunVZ, sunRadius, sunDist;
double satPara, satPerp, satDist, alpha, beta, fdummy;
int err;
meanAnomaly = tp->sunmeananomaly + ((eclTime - tp->sunepochtime) *
tp->sunmeanmotion * TWOPI);
trueAnomaly = getkep(meanAnomaly, tp->suneccentricity, &err);
satpos(tp->sunepochtime, tp->sunRAAN, tp->sunargperigee,
SUNSEMIMAJORAXIS, tp->suninclination, tp->suneccentricity, 0.0, 0.0,
eclTime, trueAnomaly, &fdummy, &fdummy, &sunX, &sunY, &sunZ, &sunRadius,
&sunVX, &sunVY, &sunVZ);
ae(sunX, sunY, sunZ, siteX, siteY, siteZ, tp->sitemat,
&tp->sunazimuth, &tp->sunelevation);
sunDist = absol(sunX, sunY, sunZ);
satDist = absol(satX, satY, satZ);
satPara = (satX*sunX + satY*sunY + satZ*sunZ) / sunDist;
satPerp = sqrt(SQR(satDist) - SQR(satPara));
alpha = asin(EARTHRADIUS / sunDist);
beta = atan(satPerp / (satPara + sunDist));
/* case 1: satellite is on side of the Earth facing the Sun */
if (satPara > 0.0) {
#ifdef ORIGINAL_CODE
if (satElevation < 0.0)
return(DAY);
if (satElevation >= 0.0 && tp->sunelevation <= TWILIGHT)
return(VISIBLE);
else
#endif
return (int) 'Y'; /* in sunlight */
}
/* case 2: satellite is on side of the Earth facing away from the Sun */
if (beta >= alpha) { /* satellite can see the Sun */
return (int) 'Y';
#ifdef ORIGINAL_CODE
if (satElevation >= 0.0 && tp->sunelevation <= TWILIGHT)
return(VISIBLE);
else
return(DAY);
#endif
}
else /* satellite is in the Earth's shadow */
return (int) 'N'; /* not in sun */
}
/*
* confss: initializes the Sun's Keplerian elements for a given
* epoch. Formulas are from "Explanatory Supplement to the
* Ephemeris, HMSO, 1974".
* It also initializes the sidereal reference when the simpler
* model is used.
*/
void confss(tp, initTime)
track_t *tp;
double initTime;
{
double T, T2, T3, omega;
int n;
if (fabs(initTime - tp->lastsuntime) <= SUNUPDATEINT)
return;
tp->lastsuntime = initTime;
/* time argument in Julian centuries since 1900.0 */
/* see Astronomical Almanac, 1983, page B6 */
/* for nutation correction see Explanatory Supplement, 1974, page 44 */
tp->sunRAAN = 0.0;
tp->sunepochtime= initTime;
tp->sidday = floor(initTime);
T = (initTime - 0.5) / 36525.0;
T2 = T*T;
T3 = T2*T;
tp->suneccentricity = 0.01675104 - 4.180e-5*T - 1.260e-7*T2;
tp->sunargperigee = (281.220844 + 1.719175*T + 4.528e-4*T2 +
3.3333e-6*T3) * CDR;
tp->sunmeanmotion = 1.0 / (365.24219878 - 6.134e-6*T);
omega = (259.183275 - 1934.142008*T + 2.078e-3*T2 +
2.22e-6*T3) * CDR;
n = (int) (omega / TWOPI);
omega -= (double) n * TWOPI;
tp->suninclination = ((23.452294 - 0.0130125*T - 1.64e-6*T2 + 5.03e-7*T3) +
2.558333e-3*cos(omega)) * CDR; /* page 41 */
#ifdef SIMPLER_MODEL
tp->sidtime = (6.646065556 + 2400.051262*T + 2.580556e-5*T2) / 24.0;
tp->sidtime -= floor(tp->sidtime);
#else
/* sia - use local apparent sidreal computed in getJulianDate */
#endif
/* for ephemeris of the Sun see Explanatory Supplement, 1974, page 98 */
tp->sunmeananomaly = (358.475833 + 35999.04975*T - 1.5e-4*T2 -
3.3333e-6*T3) * CDR;
n = (int) (tp->sunmeananomaly / TWOPI);
tp->sunmeananomaly -= (double) n * TWOPI;
tp->suntrueanomaly = getkep(tp->sunmeananomaly, tp->suneccentricity, &n);
tp->sundistance = SUNSEMIMAJORAXIS * (1.0 - SQR(tp->suneccentricity))
/ (1.0 + tp->suneccentricity * cos(tp->suntrueanomaly));
}
/*
* Theoretical free space signal attenuation. The result is *NOT* a good
* indication of the link budget. It does not take into consideration
* the effects of atmospheric composition (signal scattering and absorbtion),
* ionspheric conditions (loss due to changes in polarization, and effects
* of high energy particles, etc.), antenna apprature and gain parameters,
* etc. etc. SN101094.
*/
static float
attenuation(range, freq)
double range, freq;
{
return (float) (20.0 * log10 ( freq * ((FOURPI * range) / CVAC) ));
}
static int
ta(timeArg, tp) /* True Anomaly */
double timeArg;
track_t *tp;
{
double deltaT, fdummy, curMotion, curOrbit;
int error;
deltaT = timeArg - tp->epochday;
curMotion = tp->epochmeanmotion + deltaT * tp->decayrate;
tp->semimajoraxis = KEPLER * exp(2.0 * log(MPD / curMotion) / 3.0);
curOrbit = tp->reforbit + deltaT * curMotion;
tp->orbitnum = (long) curOrbit;
tp->meananomaly = (curOrbit - (double) tp->orbitnum) * TWOPI;
tp->trueanomaly = getkep(tp->meananomaly, tp->eccentricity, &error);
tp->orbitfrac = (int) (modf(curOrbit, &fdummy) * 100.0);
tp->orbitnum += 1;
return error;
}
void satphase(tp)
track_t *tp;
{
int n, phase;
double dphase;
dphase = (tp->meananomaly / TWOPI * tp->maxphase + tp->perigeephase);
phase = (int) dphase;
n = 0;
while (phase < 0 && ++n < 1000)
phase += (int) tp->maxphase;
n = 0;
while (phase >= (int) tp->maxphase && ++n < 1000)
phase -= (int) tp->maxphase;
tp->phase = phase;
}
/*
* satmode: gets satellite mode from phase
*/
void satmode(sp, phase, modestrp)
satellite_t *sp;
int phase;
char *modestrp;
{
mode_t *mp;
int curMode;
if (sp && (mp = sp->s_modep) ) {
for (curMode = 0; curMode < mp->nmodes; curMode++) {
if ( ( phase >= mp->mds[curMode].minPhase &&
phase < mp->mds[curMode].maxPhase) ||
( (mp->mds[curMode].minPhase > mp->mds[curMode].maxPhase) &&
( phase >= mp->mds[curMode].minPhase ||
phase < mp->mds[curMode].maxPhase) ) ) {
strncpy(modestrp, mp->mds[curMode].modeStr,2);
return;
}
}
}
strcpy(modestrp,"N/A");
}
void
predict_init(flag, sp, ptime)
int flag;
select_t *sp;
double ptime;
{
track_t *tp;
result_t *rp;
tp = sp->sl_tp;
rp = sp->sl_rp;
params_init(tp->satp, tp);
confss(tp, ptime);
siderealtime(tp, ptime);
tp->semimajoraxis = KEPLER * exp(2.0 * log(MPD / tp->epochmeanmotion) / 3.0);
tp->reforbit = (double) tp->epochorbitnum + tp->epochmeananomaly / TWOPI;
predict(flag, sp, ptime);
}
int
predict(flag, sp, ptime)
int flag;
select_t *sp;
double ptime;
{
site_t *sitep;
satellite_t *satp;
track_t *tp;
result_t *rp;
double siteX, siteY, siteZ, siteVX, siteVY, argperigee, curRaan;
double satX, satY, satZ, satVX, satVY, satVZ, satRadius, dX,dY,dZ;
double avis;
extern void satposSGP4(select_t *, double, double *, double *, double *,
double *, double *, double *, double *);
extern void SquintAngle(select_t *, double, double, double, double);
tp = sp->sl_tp;
rp = sp->sl_rp;
sitep = tp->sitep;
satp = tp->satp;
confss(tp, ptime);
siderealtime(tp, ptime);
if (fabs(tp->lastprectime - ptime) > PRECUPDATEINT) {
precess(tp->semimajoraxis, tp->eccentricity, tp->inclination,
&tp->RAANprec, &tp->perigeeprec);
tp->lastprectime = ptime;
}
sitepos(tp, sitep->c_lat, sitep->c_lng, sitep->c_alt, ptime,
&siteX, &siteY, &siteZ, &siteVX, &siteVY);
#ifdef RARELYIFEVERATALL
if (ta(ptime, tp))
return 0;
if (!satpos(tp->epochday, tp->epochRAAN, tp->epochargperigee,
tp->semimajoraxis, tp->inclination, tp->eccentricity,
tp->RAANprec, tp->perigeeprec, ptime, tp->trueanomaly, &curRaan,
&argperigee, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ))
return 0;
ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
&rp->r_elevation);
printf("MEAN X: %.3lf Y: %.3lf Z: %lf\n", satX, satY, satZ);
printf(" AZ: %.3lf EL: %.3lf\n", rp->r_azimuth, rp->r_elevation);
satposSGP4(sp, ptime, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ);
ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
&rp->r_elevation);
printf("SGP4 X: %.3lf Y: %.3lf Z: %lf\n", satX, satY, satZ);
printf(" AZ: %.3lf EL: %.3lf\n", rp->r_azimuth, rp->r_elevation);
#else
if (flag == TLEMEAN) {
if (ta(ptime, tp))
return 0;
if (!satpos(tp->epochday, tp->epochRAAN, tp->epochargperigee,
tp->semimajoraxis, tp->inclination, tp->eccentricity,
tp->RAANprec, tp->perigeeprec, ptime, tp->trueanomaly,
&curRaan, &argperigee, &satX, &satY, &satZ, &satRadius,
&satVX, &satVY, &satVZ))
return 0;
SquintAngle(sp, tp->trueanomaly, curRaan, argperigee, satRadius);
}
else if (flag == NORADSGP4)
satposSGP4(sp, ptime, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ);
else {
; /* deep-space (!) prop model. */
}
ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
&rp->r_elevation);
#endif /* _DEBUG_ */
#ifdef RANGE_USEROUTINE
satrange(siteX, siteY, siteZ, siteVX, siteVY, satX, satY, satZ,
satVX, satVY, satVZ, &rp->r_range, &rp->r_rangerate);
#else
dX = satX - siteX;
dY = satY - siteY;
dZ = satZ - siteZ;
rp->r_range = absol(dX, dY, dZ);
rp->r_rangerate = ((satVX-siteVX)*dX + (satVY-siteVY)*dY + satVZ*dZ) /
rp->r_range;
#endif /* RANGE_USEROUTINE */
ssp(tp, satX, satY, satZ, ptime, &rp->r_lat, &rp->r_long,
&rp->r_height, &avis);
rp->r_doppler = (float) ((- tp->beacon * rp->r_rangerate) / CVAC); /* in HZ */
rp->r_pathloss = attenuation(rp->r_range, tp->beacon);
satphase(tp);
satmode(satp, rp->r_phase = tp->phase, rp->r_modestr);
rp->r_insun = (char) insun(tp, satX, satY, satZ, siteX, siteY, siteZ,
ptime, rp->r_elevation);
rp->r_av = avis * CRD;
rp->r_orbitnum = tp->orbitnum;
rp->r_lat *= (float) (CRD);
rp->r_long *= (float) (CRD);
rp->r_pelev = rp->r_elevation;
rp->r_elevation *= (float) (CRD);
rp->r_azimuth *= (float)(CRD);
rp->r_time = ptime;
rp->r_rising = (rp->r_elevation >= rp->r_pelev) ? TRUE : FALSE;
return 1;
}
static void
params_init(satp, trackp)
satellite_t *satp;
track_t *trackp;
{
mode_t *mp;
trackp->satname = satp->s_name;
trackp->satnum = satp->s_satnum;
trackp->elementset = satp->s_elementset;
trackp->epochday = satp->s_epochtime;
trackp->inclination = satp->s_inclination * CDR;
trackp->epochRAAN = satp->s_raan * CDR;
trackp->eccentricity = satp->s_eccentricity;
trackp->epochargperigee = satp->s_argperigee * CDR;
trackp->epochmeananomaly = satp->s_meananomaly * CDR;
trackp->epochmeanmotion = satp->s_meanmotion;
trackp->decayrate = satp->s_decayrate;
trackp->epochorbitnum = satp->s_orbitnum;
trackp->flags = 0;
trackp->lastjdate = ZERO;
trackp->lastsuntime = ZERO;
trackp->lastprectime = ZERO;
if ((fabs(satp->s_meanmotion - 1.0) < 0.05) &&
(fabs(satp->s_inclination) < 1.0)) /* XXX FIX - sia */
satp->s_flags |= SF_GEOSTAT;
if (satp && (mp=satp->s_modep)) {
trackp->beacon = mp->beacon;
trackp->perigeephase = mp->perigeePhase;
trackp->maxphase = mp->maxPhase;
}
else {
trackp->beacon = BEACON * 1.0e6;
trackp->perigeephase = 0.9;
trackp->maxphase = MAXPHASE;
}
}
void
printMode(fp, sp, phase)
FILE *fp;
satellite_t *sp;
int phase;
{
int curMode;
mode_t *mp;
if (sp && (mp=sp->s_modep)) {
for (curMode = 0; curMode < mp->nmodes; curMode++) {
if ((phase >= mp->mds[curMode].minPhase
&& phase < mp->mds[curMode].maxPhase)
|| ((mp->mds[curMode].minPhase > mp->mds[curMode].maxPhase)
&& (phase >= mp->mds[curMode].minPhase
|| phase < mp->mds[curMode].maxPhase))) {
fprintf(fp," %2s",mp->mds[curMode].modeStr);
}
}
}
else
fprintf(fp," ");
}